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Abstract 

This  paper  presents  a  one-dimensional  steady  state  model  of  a  solid  oxide  fuel  cell  (SOFC)  to  be  used  in  an  Auxiliary  Power  Unit  (APU). 
The  fuel  cell  is  fed  a  prereformed  gas  from  an  external  autothermic  reformer.  In  addition  to  the  three  electrochemical  reactions  (reduction 
of  oxygen  at  the  cathode,  oxidation  of  hydrogen  and  carbon  monoxide  at  the  anode)  the  water-gas  shift  reaction  and  the  methane  steam 
reforming  reaction  are  taken  into  account  in  the  anode  channel.  The  model  predicts  concentrations  and  temperatures  and  uses  an  equivalent 
circuit  approach  to  describe  the  current-voltage  characteristics  of  the  cell.  The  model  equations  are  presented  and  their  implementation  into 
the  commercial  mathematical  software  FEMLAB  is  discussed.  An  application  of  this  model  is  used  to  determine  suitable  operating  parameters 
with  respect  to  optimum  performance  and  allowable  temperature. 

©  2004  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  fuel  cell  technology  will  be  an  important  part  of  the 
future  energy  supply.  Possible  applications  range  from  car 
engines  to  combined  heat  and  power  plants.  Further  on  fuel 
cells  can  be  used  to  provide  electricity  for  mobile  applica¬ 
tions.  Auxiliary  Power  Units  (APU)  are  examples  for  such  an 
application.  An  APU  based  on  a  solid  oxide  fuel  cell  (SOFC) 
is  developed  by  the  Webasto  AG.  It  is  supposed  to  supply 
trucks  with  electricity  and  it  uses  the  available  fuel  -  diesel 
or  petrol.  A  hydrogen  and  carbon  monoxide  rich  gas  is  pro¬ 
duced  in  an  external  reformer  and  fed  to  the  anode  of  the  fuel 
cell.  In  contrary  to  gases  normally  used  in  solid  oxide  fuel 
cells  this  reformated  gas  contains  less  water  and  more  than 
50%  nitrogen  due  to  the  air  utilised  in  the  reforming  process. 
At  the  cathode  pre-heated  air  is  used. 

The  electrochemical  reactions  are  described  by  an  equiva¬ 
lent  circuit  diagram.  The  use  of  equivalent  ohmic  resistances 

*  Corresponding  author.  Tel.:  +49  391  611  0350;  fax:  +49  391  611  0353. 

E-mail  address:  sundmacher@mpi-magdeburg.mpg.de 
(K.  Sundmacher). 

0378-7753/$  -  see  front  matter  ©  2004  Elsevier  B.V.  All  rights  reserved. 
doi:10.1016/j.jpowsour.2005.02.011 


for  the  mathematical  description  of  the  electrochemical  reac¬ 
tions  in  a  solid  oxide  fuel  cell  was  first  proposed  by  Achen- 
bach  [1].  This  approach  was  used  in  several  models  to  simu¬ 
late  the  steady  state  performance  of  a  SOFC  [2-4].  Recently 
there  are  endeavours  to  simulate  the  dynamic  behaviour  using 
such  an  equivalent  circuit  diagram  [5]. 

In  this  paper  a  model  of  a  solid  oxide  fuel  cell  for  an  APU 
system  is  developed  to  predict  suitable  operating  conditions. 
The  specific  composition  of  the  anode  gas  as  well  as  mate¬ 
rial  specific  limitations  are  considered.  Further  on  it  is  shown 
how  the  operating  point  of  the  fuel  cell  can  be  optimised. 
For  a  fast  and  efficient  simulation  of  the  operating  range  a 
one-dimensional  steady  state  model  is  sufficient.  The  simu¬ 
lations  are  done  with  the  commercial  mathematic  software 
FEMLAB. 


2.  Modelling 

The  model  presented  here  considers  the  smallest  repeating 
unit  of  a  co-flow  or  counter-flow  solid  oxide  fuel  cell  stack. 
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Nomenclature 


Latin  letters 


c 

concentration  (mol  m-3) 

cv 

molar  heat  capacity  (J  mol-1  K- 

Ea 

activation  energy  (J  mol-1) 

F 

Faraday’s  constant  (A  s  mol-1) 

g 

molar  flux  density  (molm-2  s-1 

Ar  g 

enthalpy  of  reaction  (J  mol-1) 

h 

height  (m) 

I 

current  (A) 

i 

current  density  (Am-2) 

k 

reaction  rate  constant  (molm-3) 

K 

equilibrium  constant 

n 

number  of  transferred  electrons 

Ar  h 

heat  of  reaction  (J  mol-1) 

P 

pressure  (Pa) 

91 

heat  flux  density  (J  m-2  s-1) 

R 

gas  constant  (J  mol-1  K-1) 

R 

area-related  resistance  (Q  m2) 

r 

volumetric  reaction  rate  (mol  m- 

T 

temperature  (K) 

U 

voltage  (V) 

f^Nernst 

Nernst  voltage  (V) 

U 

velocity  (ms-1) 

X 

molar  fraction 

z 

coordinate  (m) 

Greek  letters 

a 

heat-transfer  coefficient  (W  m-2 

X 

heat  conductivity  (W  m-1  K-1) 

V 

stoichiometric  coefficient 

Subscripts 

a 

anode 

avg 

average 

c 

cathode 

el 

electric 

eq 

equilibrium 

FC 

fuel  cell 

i 

index  for  the  component 

in 

inlet 

3 

index  for  the  reaction 

pol 

polarisation 

ref 

methane  reforming 

t 

total 

wgsr 

water-gas  shift  reaction 

Components 

ch4 

methane 

CO 

carbon  monoxide 

co2 

carbon  dioxide 

h2 

hydrogen 

—  l 


H2O  water 

N2  nitrogen 

O2  oxygen 


It  consists  of  a  single  anode  and  cathode  channel  (Fig.  1). 
Thereby  the  complexity  of  the  problem  is  reduced  and  the  nu¬ 
merical  computation  time  is  short  enough  to  perform  repeated 
simulation  in  an  acceptable  time  frame,  making  optimisation 
possible.  Furthermore,  the  model  complexity  is  reduced  to 
only  one  spatial  dimension  by  assuming  plug  flow  in  both 
channels,  that  is  setting  all  gradients  along  the  x-  and  y- 
coordinate  to  zero. 

2.7.  Model  assumptions 

The  smallest  repeating  unit  of  a  SOFC  stack  consists  of 
three  compartments:  the  cathode  and  the  anode  channel  as 
well  as  the  solid  phase.  The  solid  phase  comprises  all  solid 
parts  of  the  model  (the  bipolar  plate,  the  cathode  and  anode 
electrode  and  the  electrolyte).  The  structure  of  the  model  is 
shown  in  Fig.  2. 


concentration  temperature  reactions 

velocity  heat  and  mass  transfer 


(Tc) 

cathode  channel 

C0,.c  Tc  uc  |q«c  |S‘-C 

(Tc) 
.... _ 

(Ts) 

solid  phase  -p  ,  .... 

1  s  electrochemical  reactions 

1  1 

(Ts) 

<Ta> 

CH,.a  CHAa  Ta  Ua  \Q  l„ 

r*  n  “a.  a  ©i.a 

a  '-fO  a 

c  :  ‘  methane  reforming 

CH4  ,a 

,  water  gas  shift  reaction 

anode  channel 

_ 

O’.) 

z=o  7.  =  1, 


co-flow 

counter-flow 
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Fig.  2.  Schematic  one-dimensional  representation  of  the  SOFC  model. 
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The  temperature  is  calculated  for  each  compartment  of 
the  model.  Further  on  the  velocity  and  the  concentration  of 
the  components  are  considered  for  the  cathode  and  the  anode 
channel.  The  gas  phases  and  the  solid  phase  are  intercon¬ 
nected  via  exchange  of  mass  (molar  flux  density,  gi)  and  heat 
(heat  flux  density,  qa). 

The  electrochemical  reactions  take  place  inside  the  porous 
electrodes,  which  are  considered  to  be  a  part  of  the  solid 
phase.  At  the  cathode  oxygen  is  reduced 

02  +  4e“  — 202-  (1) 

whereas  at  the  anode  hydrogen  and  carbon  monoxide  are 
oxidised 

H2  +  02“  -*  H20  +  2e“  (2) 

CO  +  O2-  -*  C02  +  2  e“  (3) 

The  reactions  mentioned  in  Eqs.  (l)-(3)  result  in  the  fol¬ 
lowing  two  overall  reactions: 

H2  +  ±02+  -*  H20  (4) 

CO  +  \  02+  -*  C02  (5) 

The  fuel  is  converted  to  a  hydrogen  and  carbon  monoxide 
rich  gas  in  an  external  reformer,  which  is  not  part  of  the 
model.  The  resulting  gas  contains  a  small  amount  of  methane, 
which  undergoes  the  strongly  endothermic  methane  steam 
reforming  reaction  in  the  anode  channel: 

CH4  +  H20  ^  CO  +  3  H2  (6) 


neglected  in  the  gas  phases.  In  the  solid  phase  heat  con¬ 
duction  along  the  spatial  coordinate  is  considered. 

•  The  methane  steam  reforming  reaction  and  the  water- 
gas  shift  reaction  are  considered  as  reversible,  quasi- 
homogenous  reactions  which  are  at  the  chemical  equilib¬ 
rium.  The  enthalpy  of  this  reactions  is  taken  into  account 
in  the  energy  balance  of  the  anode  channel. 

•  The  electrical  conductivity  of  the  bipolar  plate  is  suffi¬ 
ciently  high;  therefore  the  ohmic  resistance  can  be  ne¬ 
glected.  As  a  result  the  cell  voltage  is  constant  along  the 
spatial  cell  coordinate. 

•  There  are  no  concentration  polarisations  due  to  internal 
mass  transport.  In  reality  these  effects  are  only  observable 
at  low  concentrations  of  an  educt  component  or  at  high 
current  densities. 

2.2.  Model  equations 

Mass  balances  in  both  gas  phases  describe  the  composi¬ 
tion  and  velocity  in  the  anode  and  cathode  channel.  Enthalpy 
balances  in  temperature  form  yield  the  temperature  in  the 
gas  phases  and  in  the  solid  phase.  Kinetic  expressions  are 
required  for  the  heat  exchange  between  the  phases  and  for 
the  reaction  rates.  The  thermodynamic  expressions  from  the 
National  Institute  of  Science  and  Technology  [6]  are  used  to 
describe  the  temperature  dependence  of  the  gas  properties. 
The  coefficient  equations  for  the  heat  capacity,  the  enthalpy 
and  the  Gibbs  enthalpy  are  applied  in  the  model. 


Further  on  the  slightly  exothermic  water-gas  shift  reaction  is 
considered  in  the  anode  gas  phase: 

CO  +  H20  ^  C02  +  H2  (7) 

Both  reactions  are  included  in  the  homogenous  model  of  the 
gas  in  the  anode  channel. 

Due  to  the  design  of  the  bipolar  plate  there  is  an  inactive 
region  before  and  after  the  active  cell  area  (see  Fig.  2).  Within 
this  region  only  the  heat  exchange  between  the  anode  gas, 
the  solid  and  the  cathode  gas  takes  place.  Neither  the  gas 
phase  reactions  in  the  anode  channel  nor  the  electrochemical 
reactions  are  considered  here. 

Further  on  the  following  assumptions  are  applied  for  the 
mathematical  derivation  of  the  model: 

•  Steady  state  operation  is  considered. 

•  The  ideal  gas  law  is  applied  for  all  gas  components. 

•  The  pressure  in  the  model  is  considered  as  constant  and 
equal  to  the  standard  pressure  p°  =  101.325  kPa.  This 
condition  is  approximately  fulfilled  in  the  real  world  by  the 
respective  system  design,  which  considers  ambient  pres¬ 
sure  and  low  pressure  drop  along  the  channel. 

•  Plug-flow  is  assumed  within  the  cathode  and  the  anode 
channel.  Thereby  the  gas  temperature,  the  gas  velocity 
and  the  composition  only  change  along  the  z-coordinate. 

•  Only  convective  mass  and  heat  transfer  in  z-direction  are 
considered;  diffusion  and  heat  conduction  phenomena  are 


2.2.7.  Balance  equations 

The  balance  equations  for  both  gas  channels  are  similar. 
Therefore,  only  the  equations  for  the  anode  channel  are  de¬ 
scribed  below.  The  equations  for  the  cathode  channel  can 
be  obtained  analogously.  The  only  significant  difference  be¬ 
tween  them  is  that  no  gas  phase  reaction  occurs  in  the  cathode 
channel. 

The  differential  equation  for  the  concentration  of  the  com¬ 
ponents  in  the  gas  phases  is  derived  from  the  steady  state 
continuous  mass  balance. 


0  =  - 


d 

3  z 


(ci,  a  ua)  +  7^ - F  Vi, 

d a  j 


j  rj 


Besides  the  convection  term  the  mass  transfer  from  and  to  the 
solid  phase  and  the  quasi-homogenous  reactions  in  the  gas 
phase  are  covered  by  this  equation. 

The  following  components  i  and  reactions  j  are  considered 
for  the  anode  channel  (index  a): 


i  =  H2,  CO,  C02,  CH4,  H20,  j  =  ref,  wgsr 

For  the  cathode  channel  (index  c)  only  the  concentration  of 
oxygen  is  calculated.  There  are  no  reactions  in  this  channel. 


i  =  02, 


7  —  0 
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The  concentration  of  nitrogen  is  obtained  from  the  summa¬ 
tion  condition  ct  =  JA  Ci  using  the  ideal  gas  law  ct  =  p/(RT) 
for  the  total  concentration. 


cn2  = 


The  differential  equation  for  the  velocity  of  the  gas  is  de¬ 
rived  from  the  overall  mass  balance: 


0  =  - 


3 

3  z 


(ct,aMa)  +  J2  + 

na 


VJ  ri 


(10) 


Insertion  of  the  ideal  gas  law  into  Eq.  (10)  yields 


0  =  - 


p  3/1 


RdzKZ 


u 


na 


a 


a 


(ID 


To  solve  Eq.  (11)  an  additional  differential  equation  for 
the  temperature  ra  must  be  given.  The  needed  equation  is 
obtained  from  the  steady  state  enthalpy  balance: 


0  =  - 


3T 


uact,a  cp 


a 


,a 


3  Z 


+  y~l(— Ar  hj)fj  + 
j 


Ola, a 


(12) 


with  the  average  heat  capacity  c p  =  JA(cp,z<7/ct)- 

The  first  term  on  the  right  hand  side  of  Eq.  (12)  shows 
the  convective  enthalpy  transport  whereas  the  second  term 
shows  the  temperature  change  as  a  result  of  the  mass  transfer 
between  the  solid  phase  and  the  gas  phase.  In  Eq.  (12),  only 
the  gas  flow  from  the  solid  phase  into  the  gas  channel  is  taken 
into  account  (see  Fig.  2)  [7]: 


Ta,  ifg/,a<0 
Ts,  ifgz',a  >  0 


(13) 


Further  on,  the  reaction  heat  of  the  gas  phase  reactions  in 
the  anode  channel  and  the  heat  transfer  between  the  two  com¬ 
partments  are  included.  A  linear  approach  is  used  to  describe 
the  heat  exchange: 

Ola, a  —  Oia{Ts  ~  7a).  (14) 


The  value  of  the  heat-transfer  coefficient  aa  is  calculated 
for  a  reference  velocity  and  composition  in  the  anticipated 
temperature  range  [8].  A  linear  function  is  estimated  from 
these  values  to  describe  the  temperature  dependence  of  the 
heat- transfer  coefficient  in  the  model. 

A  differential  equation  for  the  temperature  of  the  solid 
phase  can  be  obtained  from  the  steady  state  enthalpy  balance. 
Heat  conduction  according  to  Fourier’s  law,  thermal  drift  due 
to  mass  transfer  from  and  to  the  gas  phase,  convective  heat 
exchange  and  a  heat  source  due  to  the  electrochemical  reac¬ 
tions  balance  each  other: 

3  2rs 

0  =  hsXs — tt- 
dz2 

T  ^  ^  6p,i  §i, c(7c/s  —  7s)  —  Ola, a  ~  Ola,c  T  Ols  (13) 


Table  1 

Boundary  conditions  applied  to  the  differential  Eqs.  (8),  (11),  (12)  and  (15) 
for  co-flow  and  counter-flow  configuration 


Left  hand  side  (z  =  0) 

Variable 

Right  hand  side  (z  =  L) 

Cathode  channel  (co-flow) 

C02,c  —  C02,c,in 

U  c  —  Wcin 

C02,c 

uc 

dC(Vc  _  n 

dz  _  U 

d."c  =  0 
dz 

Tc  —  Tc,  in 

Tc 

a 

*  o* 

II 

o 

Cathode  channel  (counter-flow) 

dc02,c  _  n 

dz  “  U 

C02,c 

C02,c  —  C02,c,in 

d«c  _  n 

dz  ~  U 

Uc 

Uc  —  Uc,m 

o 

II 

£1- 

Tc 

Tc  =  Tcjn 

Solid  phase 

f  =0 

Ts 

II 

o 

Anode  channel 

Ta  —  Ta[n 

—  ^a,in 

Ta 

ua 

^  =  0 

dwa  _  n 

dz  “  U 

cH2,a  =  CH2,a,in 

cH2,a 

dcH2,a  n 

difz  —  J 

CCO,a  =  CC0,a,in 

CCO,a 

dcCO,a  n 

dz  _  U 

cC02,a  =  CC02,a,in 

cC02,a 

dcco2,a  p. 

dz  =0 

CH20,a  —  QrbCyadn 

cH20,a 

dfHoO.a  p. 

dz  =0 

CCH4,a  =  CCH4,a,in 

CCH4,a 

dccH4,a  n 

dz  “  U 

FEMLAB  requires  two  boundary  conditions  for  each  ODE. 


The  enthalpy  of  both  overall  reactions  (see  Eqs.  (4)  and 
(5))  are  taken  into  account.  A  part  of  this  energy  is  converted 
into  electric  power.  The  remaining  energy  is  released  as  heat: 

qs  =  (— AR/zel,H2)£H2  +  (-  AR/zel,Co)gCO  “  ^FC  *el-  (16) 

In  Eq.  (16),  —  Ar/z^  h?  and  —  Ar/z^co  characterise  the 
reaction  heat  of  the  oxidation  of  hydrogen  and  carbon  monox¬ 
ide,  respectively.  gn2  and  geo  describe  the  molar  flux  densi¬ 
ties  between  the  gas  channel  and  the  solid  and  can  be  consid¬ 
ered  as  a  representation  of  the  reaction  rate  per  unit  electrolyt 
area. 

For  each  of  the  presented  differential  equations,  two 
boundary  conditions  are  needed  due  to  the  fact  that  the  sim¬ 
ulation  software  FEMLAB  uses  parabolic  differential  equa¬ 
tions  of  second  order  with  respect  to  the  spatial  coordinate 
(Table  1).  The  inlet  properties  of  the  anode  gas  are  assigned 
to  the  left  hand  side  of  the  model  (z  =  0).  According  to  the 
flow  direction  the  inlet  velocity,  temperature  and  concentra¬ 
tion  for  the  cathode  channel  are  specified  at  the  left  hand 
side  (co-flow:  z  =  0)  or  at  the  right  hand  side  (counter-flow: 
z  =  L).  For  the  second  boundary  a  zero  gradient  is  assumed. 
Perfect  insulation  is  assumed  for  both  boundaries  of  the  solid 
phase. 

2.2.2.  Kinetics  of  the  anode  gas  phase  reactions 

In  the  anode  channel  of  the  model  the  methane  steam  re¬ 
forming  reaction  and  the  water-gas  shift  reaction  take  place. 
Different  expressions  for  the  kinetics  of  both  reactions  can 
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be  found  in  the  literature  (for  example  [2,3,9]).  The  kinetics 
depend  on  the  exact  composition  of  the  electrode  or  catalyst 
material.  They  are  not  assignable  to  our  simulations  with¬ 
out  great  effort.  Therefore,  simple  reaction  mechanisms  are 
used. 

In  this  model  chemical  equilibrium  of  the  gas  phase  reac¬ 
tions  is  assumed.  According  to  Eq.  (6),  the  equilibrium  of  the 
methane  steam  reforming  reaction  is  given  by 


*CO,a-*H2,a 

•XCH4,a*H20,a 


(17) 


Analogously,  the  equilibrium  of  the  water-gas  shift  reaction 
(Eq.  (7))  can  be  written  as: 


^  _  *C02,a*H2,a 

Awgsr,eq  — 

*C0,a*H20,a 


(18) 


with  the  molar  fraction  jcz-  =  c;  /ct  and  the  equilibrium  con¬ 
stant 


Ri,c  q  —  exP 


(19) 


where  Ar g°  is  the  standard  Gibb’s  enthalpy  of  the  reaction. 

Eqs.  (17)  and  (18)  are  algebraic  equations,  thus  the  result¬ 
ing  equation  system  is  a  differential  algebraic  system  (DAE). 
The  implementation  of  the  equilibrium  condition  in  algebraic 
form  leads  to  unacceptable  numerical  complication.  There¬ 
fore,  we  substitute  these  equations  with  finite  reaction  kinet¬ 
ics  which  fulfil  the  thermodynamic  equilibrium  conditions. 
To  achieve  gas  compositions  close  to  equilibrium  at  any  point 
within  the  active  cell  area,  the  reaction  rate  constants  are  arbi¬ 
trarily  set  to  high  values.  Simple  power-law  reaction  kinetics 
with  an  Arrhenius  term  are  used  for  both  reactions: 


7] ref  —  ^ref,+,oo 


E  A, ref  \ 

RTz  J 


x 


•*CH4,a*H20,a 


*CO,a* 


3 

H2,a 


(20) 


r wgsr  —  ^wgsr,+,oo  exp  ( 


Ea, 


wgsr 


RT 


a 


X 


*C0,aAH20,a 


*C02,aAH2,a 


K 


wgsr,eq 


(21) 


2.2.3.  Electrochemical  reactions 

Two  paths  are  considered  for  the  electrochemical  reac¬ 
tions.  Either  hydrogen  or  carbon  monoxide  can  be  electro- 
chemically  oxidised  at  the  anode  electrode.  The  kinetics  of 
this  reactions  are  described  by  using  the  approach  of  Achen- 
bach  [1].  Equivalent  electric  resistances  are  defined  for  the 
overpotentials  at  the  cathode  and  at  the  anode.  Furthermore, 
it  is  assumed  that  the  resistance  is  independent  of  the  local 
current  density.  The  equivalent  circuit  diagram  given  in  Fig.  3 


Fig.  3.  Equivalent  circuit  diagram. 


illustrates  this  approach. 


7?pol,c,02 


4  F 
~RT 


ko2  (*02,c)m°2 


exp 


E  A,pol,Q2 

RTS 


(22) 


7?pol,a,H2  — 


IF 


(23) 


2  F  /  \nico  f  ^A,pol,Co\ 
~kco  (*co,a)  exp  ( - fjr-  J 

(24) 

The  constants  used  to  calculate  the  polarisation  resistances 
(Eqs.  (22)-(24))  are  listed  in  Table  2. 

The  total  cell  current  /pc  is  distributed  along  the  z- 
coordinate  as  the  local  cell  current  density  ze b  which  is 
divided  into  /ei,H2  and  zei,co  according  to  the  two  parallel 


7^pol,a,CO  — 


Table  2 


Constants  used  for  the  cathode  and  anode  polarisation  resistance  [1] 


k  (A/m2) 

£A,poi  (J/molK) 

m 

^pol,c,02 

14.9  x  109 

160  x  103 

0.25 

^pol,a,H2 

0.213  x  109 

110  x  103 

0.25 

^pol,a,CO 

0.298  x  109 

110  x  103 

0.25 
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reaction  paths: 


with  /avg  as  the  average  current  density  and  Acen  as  the  active 
cell  area.  The  electric  resistances  R  and  the  Nernst  voltages 
^Nernst  are  also  local  values  whereas  the  cell  voltage  U^c  is 
constant  along  the  channel. 

A  mathematical  description  of  the  equivalent  circuit  dia¬ 
gram  is  obtained  from  Kirchhoff’s  laws: 


hi  =  4i,h2  +  *el,CO  (26) 

Un2  =  UCo  (27) 

£4t2  =  ^Nernst,H2  “  *el,H2  '  ^pol,a,H2  (28) 

UcO  =  ^Nernst, CO  “  hi, CO  '  ^pol,a,CO  (29) 

Ufc  =  Uh2  -  (Rpoic  +  Rn)hi  (30) 

Rq  =  Rq,  c  +  RqI  +  Rq,sl  (31) 


whereas  the  ohmic  resistance  of  the  anode  electrode  (/to,a) 
and  the  cathode  electrode  (Rq,c)  as  well  as  the  resistance  for 
the  oxygen  ions  in  the  electrolyte  (RQ\)  are  calculated  using 
the  equations  given  by  Achenbach  [1]. 

The  Nernst  voltage  required  in  the  Eqs.  (28)  and  (29)  is 
given  by: 


f^Nernst,H2 


—  ARg°H2  RT^  ln  / V*o2,c*H2,a\ 
2 F  2F  n  V  *H2o,a  ) 


and 


f^Nernst,CO 


-ARg°co  RT^  ( V*02,c*CO,a\ 

IF  2 F  \  vco2,a  )  ' 


(32) 


(33) 


The  connection  between  the  current  density  and  the  molar 
flux  density  is  given  by  Faraday’s  law.  It  is  applied  for  the 
oxidation  of  hydrogen  and  the  oxidation  of  carbon  monoxide 
at  the  anode  and  analogously  for  the  reduction  of  oxygen  at 
the  cathode. 


*el,H2  =  ^H2gH2,a^  (34) 

*’ei,co  =  ncogco,&F  (35) 

hi  =  no2  8o2,cF  (36) 

In  the  equations  above  n  is  the  number  of  transferred  electrons 
in  the  reduction  and  oxidation  process,  respectively.  The  mo¬ 
lar  flux  density  of  the  mass  exchange  flux  between  the  gas 
channel  and  the  solid  phase  is  represented  by  g  (see  Fig.  2). 
F  stands  for  Faraday’s  constant. 


2.3.  Definition  of  the  operating  point 

An  additional  input  parameter  is  required  to  define  the 
operating  point  in  the  current- voltage  plain  of  the  fuel  cell. 
On  the  one  hand  it  is  possible  to  specify  a  certain  cell  voltage 


(potentiostatic  mode),  on  the  other  hand  the  cell  current  can 
be  given  (galvanostatic  mode). 

For  the  potentiostatic  mode  the  cell  current  is  easily  cal¬ 
culated  from  the  charge  balance  as  the  integral  of  the  current 
density  over  the  active  cell  area  (Eq.  (25)).  An  integral  con¬ 
dition  has  to  be  solved  for  the  galvanostatic  mode: 

F\(U?c)  =  7fc  -  7*(f/pc)  =>  0  (37) 

where  /pc  is  the  given  cell  current  and  /*  is  defined  as 

7*  =  [  hi(UFc,xit&(z),XifC(z),T$(z))bdz  (38) 

Jz=h 

This  means  that  £/fc  has  to  be  varied  in  such  a  manner  that 
Eq.  (37)  is  fulfilled. 


3.  Implementation  in  FEMLAB 

The  commercial  modelling  environment  FEMFAB  2.3 
is  used  to  solve  the  set  of  governing  equations.  The 
program  provides  a  number  of  predefined  partial  differ¬ 
ential  equations  (PDEs)  from  several  areas  of  science 
and  engineering  (referred  to  as  application  modes).  Ad¬ 
ditional  model  specific  PDEs  can  be  defined.  All  equa¬ 
tions  are  solved  simultaneously  by  applying  the  finite 
element  method  (FEM).  A  specific  documentation  of 
all  abilities  of  FEMFAB  can  be  found  in  the  manuals 
[10,11]. 

A  one-dimensional  geometry  consisting  of  three  subdo¬ 
mains  is  created  by  using  the  graphical  user  interface  (GUI). 
The  middle  subdomain  is  the  active  cell  area  with  inactive 
regions  on  the  left  and  on  the  right. 

The  model  uses  differential  equations  to  calculate  the 
composition  and  velocity  of  the  gas  in  both  channels  as 
well  as  the  temperature  in  all  three  compartments.  A  con¬ 
vection  and  diffusion  application  mode  is  used  for  the 
partial  mass  balance  in  the  gas  phases  (Eq.  (8))  whereas 
the  temperatures  in  the  gas  channels  (Eq.  (12))  are  de¬ 
scribed  in  a  convection  and  conduction  mode.  In  the  solid 
phase  the  enthalpy  balance  (Eq.  (15))  is  implemented  in 
a  heat  transfer  application  mode.  Given  that  the  equation 
used  for  the  velocity  (Eq.  (11))  does  not  exist  as  a  pre¬ 
defined  application  mode,  a  pde  general  form  is  applied. 
The  boundary  conditions  for  all  equations  are  listed  in 
Table  1. 

The  integral  condition  (Eqs.  (37)  and  (38))  to  set  the  cell 
current  in  the  galvanostatic  mode  is  implemented  using  a 
weak  boundary  and  a  weak  subdomain  application  mode. 
Due  to  this,  the  weak  form  of  the  PDEs  has  to  be  used  to 
solve  the  problem. 

Furthermore,  all  equations  regarding  the  calculation  of 
interim  values  are  supplied  as  expressions.  Parameters  and 
other  constant  values  are  entered  as  constants. 

The  results  are  calculated  using  the  stationary  non-linear 
solver.  Due  to  the  strong  non-linearity  of  the  problem,  it  is 
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Table  3 

Definition  of  the  reference  case 


Description 

Variable 

Value 

Electrical  characteristics 

Current  (A) 

he 

0.45 

Gas  properties  at  the  cathode  inlet 

Temperature  (K) 

Tc,  in 

923 

Velocity  (ms-1) 

He, in 

0.4 

Compound 

Oxygen 

x02,c,'m 

0.21 

Nitrogen 

•^N2,c,in 

0.79 

Gas  properties  at  the  anode  inlet 

Temperature  (K) 

Fa,  in 

1123 

Velocity  (ms-1) 

^a,in 

0.5 

Compound 

Hydrogen 

•XH2,a,in 

0.22 

Carbon  monoxide 

•^CO,a,in 

0.20 

Carbon  dioxide 

•XC02,a,in 

0.02 

Water 

•^t^O.a.in 

0.03 

Methane 

•^CH4,a,in 

0.001 

Nitrogen 

•^N2,a,in 

0.529 

hard  to  find  an  initial  guess  such  that  the  solution  converges. 
Therefore  the  parametric  solver  is  applied.  First,  the  steady 
state  is  calculated  with  a  current  of  /pc  =  0.01  A  as  a  param¬ 
eter.  Then,  the  cell  current  is  increased  stepwise.  For  each 
new  value  the  solution  of  the  old  parameter  is  used  as  an  ini¬ 
tial  guess.  The  result  is  the  current  voltage  characteristic  of 
the  cell. 


4.  Results 

The  model  is  used  to  estimate  the  steady  state  of  the  SOFC 
stack.  Calculations  are  done  for  the  co-flow  and  the  counter¬ 
flow  configuration  and  a  first  analysis  of  the  operating  range 
is  performed. 

4.1.  Operation  conditions  at  reference  case 

The  values  used  in  the  reference  case  are  shown  in  Table  3. 
The  galvanostatic  operation  mode  is  applied  for  the  simula¬ 
tions.  The  current  of  /pc  =  0.45  A  is  equivalent  to  an  average 
current  density  of  /avg  =  250mA/cm2. 

At  the  cathode  inlet  pre-heated  air  at  a  temperature  of 
923  K  enters  the  fuel  cell  whereas  at  the  anode  inlet  gas  from 
the  external  reformer  is  used.  The  composition  approximates 
the  equilibrium  molar  fractions  after  the  partial  oxidation  of 
diesel  at  1 123  K  with  an  air  to  fuel  ratio  of  0.4. 

4.2.  Co -flow 

The  flow  direction  in  the  cathode  and  the  anode  channel 
are  identical  in  the  co-flow  case.  The  temperature  distribution 
in  both  channels  and  in  the  solid  phase  is  shown  in  Fig.  4. 
Within  the  inactive  regions,  the  temperatures  of  the  cathode 
air  and  the  anode  gas  approach  the  solid  temperature.  It  re- 


Fig.  4.  Temperature  distribution  (co-flow)  at  reference  conditions. 


suits  in  an  average  temperature  of  approximately  1035  K  at 
the  boundary  between  the  inactive  region  and  the  active  area. 
After  that  the  temperature  increases  as  a  result  of  the  heat 
which  is  released  by  the  electrochemical  reaction.  Since  the 
solid  phase  is  heated  by  the  electrochemical  reaction,  its  tem¬ 
perature  is  highest,  while  the  temperatures  of  the  cathode  air 
and  the  anode  gas  are  up  to  10  K  lower.  A  temperature  of 
approximately  1210  K  is  reached  at  the  end  of  the  active  cell 
area. 

Fig.  5  shows  the  corresponding  current  density.  The  graph 
follows  the  temperature  distribution  in  the  active  cell  area. 
The  electric  resistance  decreases  with  increasing  tempera¬ 
ture.  Due  to  Ohm’s  law  and  the  constant  cell  voltage,  the 
current  density  rises  accordingly. 

The  composition  of  the  gas  mixture  in  the  cathode  channel 
is  shown  in  Fig.  6.  The  consumption  rate  of  oxygen  increases 
with  rising  current  density.  Therefore,  the  concentration  of 
O2  in  the  cathode  channel  decreases  accordingly.  Under  the 
conditions  defined  in  the  reference  case,  the  molar  fraction 
of  oxygen  in  the  gas  decreases  from  21  to  16.5%.  The  curve 


Fig.  5.  Current  density  distribution  (co-flow)  at  reference  conditions. 
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Fig.  6.  Composition  of  the  cathode  gas  (co-flow)  at  reference  conditions. 

corresponds  to  the  temperature  and  current  density  distribu¬ 
tion. 

In  addition  to  the  electrochemical  reactions  the  water-gas 
shift  reaction  and  the  methane  steam  reforming  reaction  are 
regarded  in  the  anode  channel  (Fig.  7). 

The  concentration  of  hydrogen  is  almost  constant  in  the 
first  part  of  the  active  cell  area.  The  hydrogen  oxidised  by 
the  electrochemical  reaction  is  balanced  by  the  water-gas 
shift  reaction.  Since  the  anode  gas  is  cooled  down  within  the 
inactive  region,  the  water-gas  shift  reaction  is  no  longer  at 
the  chemical  equilibrium.  Thus  carbon  monoxide  reacts  with 
water  to  carbon  dioxide  and  hydrogen.  The  change  of  the 
composition  by  the  progressive  oxidation  of  H2  and  CO  as 
well  as  the  increase  of  the  temperature  leads  to  a  crossing  of 
the  equilibrium  of  the  water-gas  shift  reaction.  Afterwards 
the  reverse  water  gas  shift  reaction  outbalances  this  reaction. 
An  indicator  is  the  increasing  water  concentration  and  the 
relatively  stable  molar  fraction  of  carbon  monoxide  in  the 
second  half  part  of  the  active  cell  area. 


Fig.  7.  Composition  of  the  anode  gas  (co-flow)  at  reference  conditions. 


Table  4 

Simulation  results  for  the  reference  case  (co-flow) 


Simulated  state 

Variable 

Value 

Average  current  density  (mAcm2) 

*el,avg 

250.0 

Voltage  (V) 

Uvc 

0.67 

Average  power  density  (mW  cm2) 

Pel,  avg 

167.5 

Average  solid  temperature  (K) 

Ts,avg 

1116 

Maximum  solid  temperature  (K) 

Fs,max 

1213 

A  summary  of  characteristic  process  values  for  the  co-flow 
operating  mode  can  be  found  in  Table  4. 

4.3.  Counter-flow 

The  cathode  air  inlet  and  fuel  inlet  are  on  opposite  sides  for 
the  counter-flow  case.  The  resulting  temperature  distribution 
is  displayed  in  Fig.  8. 

The  gases  at  the  outlets  are  cooled  down  by  the  gas  at 
the  corresponding  inlet.  The  heat  produced  by  the  electro¬ 
chemical  reaction  can  not  be  discharged  and  the  temperature 
adjusts  to  the  displayed  distribution  for  the  steady  state.  The 
maximum  temperature  rs  max  ~  1400  K  can  be  found  in  the 
middle  of  the  active  cell  area. 

The  maximum  temperatures  are  very  high  compared  to  the 
co-flow  case.  Accordingly  this  results  in  a  higher  temperature 
stress  in  the  materials.  Due  to  the  allowed  maximal  temper¬ 
ature  of  1250  K  the  counter-flow  mode  is  not  recommended 
for  SOFC  operating  and  is  therefore  not  considered  for  the 
following  analyses  [5]. 

4.4.  Isothermal  U-I  characteristics 

The  operational  behaviour  of  the  solid  oxide  fuel  cell  stack 
is  described  by  its  current-voltage  characteristics.  The  com¬ 
putation  of  these  characteristics  can  be  done  on  the  one  hand 
by  the  specification  of  the  cell  voltage  (potentiostatic)  or  on 
the  other  hand  by  the  specification  of  the  cell  current  (gal- 


Fig.  8.  Temperature  distribution  (counter-flow)  at  reference  conditions. 


M.  Pfajferodt  et  al.  /  Journal  of  Power  Sources  149  (2005)  53-62 


61 


^Ue|1(T  =  1000K) 
Ucell  (T  =  1 100K) 
Uceli  (T  =  1200K) 
—  Ue|,  (T  =  1300K) 
o  -Pcell  (T  =  1 OOOK) 
°  PCen  (T  =  1 100K) 
*  pcell  (T  =  1200K) 
P.^„  (T-13QOK) 


Fig.  9.  Isotherm  U-I  characteristic. 


vanostatic).  Both  variants  result  in  the  same  current- voltage 
characteristics. 

The  current-voltage  characteristic  of  a  solid  oxide  fuel 
cell  is  measured  at  constant  cell  temperature.  To  make  these 
simulations  comparable  to  these  experiments,  isothermal 
conditions  are  used.  Fig.  9  shows  the  U-I  characteristic 
for  a  constant  cell  temperature  of  1000,  1100,  1200  and 
1300  K.  The  galvanostatic  mode  is  used  to  calculate  the 
curves. 

The  electric  resistance  only  depends  on  the  temperature 
and  the  concentration  of  the  educt  gases.  As  an  isothermal 
model  is  considered  and  the  influence  of  the  educt  concen¬ 
tration  shows  only  at  high  current  densities,  the  voltage  loss 
rises  linearly  according  to  Ohm’s  law.  Furthermore  the  differ¬ 
ence  between  the  curves  diminishes  with  increasing  temper¬ 
ature  and  the  maximum  of  the  power  density  shifts  to  higher 
current  densities. 

Fig.  10  shows  the  voltage  losses  caused  by  each  electri¬ 
cal  and  electrochemical  element  at  a  temperate  of  1200  K. 
The  losses  due  to  the  polarisation  resistance  at  the  cathode 
and  the  ion  transport  through  the  electrolyte  are  of  similar 
order  of  magnitude.  At  the  anode  the  oxidation  of  hydrogen 
and  the  oxidation  of  carbon  monoxide  have  to  be  considered 
together.  The  overall  voltage  loss  at  the  anode  can  be  esti¬ 
mated  with  a  parallel  connection  of  the  corresponding  polar¬ 
isation  resistances.  This  results  in  a  voltage  loss  of  about  ^ 
of  the  cathodic  loss.  The  effect  of  the  ohmic  resistances  of 
the  cathode  and  the  anode  material  on  the  cell  voltage  can  be 
neglected. 

Further  investigation  of  the  polarisation  resistances  as  well 
as  an  analysis  of  sensitivity  to  the  thickness  of  the  individual 
cell  components  are  done  by  Chen  et  al.  [12].  They  compare 
activation  and  concentration  polarisation  for  electrolyte  sup¬ 
ported  and  anode  supported  cells.  The  results  for  electrolyte 
supported  cells  correspond  qualitatively  to  the  results  shown 
in  Fig.  10. 

4.5.  Operating  range 

The  solid  oxide  fuel  cell  used  for  the  APU  system  repre¬ 
sents  a  complex  chemical  and  physical  system.  Due  to  the 
material  stresses  there  are  limitations  to  the  operating  range, 


i.e.  there  are  restrictions  for  the  system  variables.  The  relevant 
variables  are  the  concentrations,  the  velocities,  the  current 
density  and  the  temperature  distributions. 

The  presented  simulation  describes  the  steady  state  of  the 
solid  oxide  fuel  cell.  The  steady  state  of  the  system  depends 
on  the  value  of  the  input  parameters,  which  can  be  divided 
into  three  groups  according  to  their  impact: 

•  The  fuel  at  the  anode  inlet: 

The  temperature  (ra  in),  the  velocity  (wa, in)  as  well  as 
the  composition  of  the  anode  gas  (x^ in)  have  to  be  de¬ 
fined  for  the  simulation.  The  needed  hydrogen  and  carbon 
monoxide  rich  gas  is  produced  by  the  partial  oxidation  of 
diesel  in  the  reformer.  For  the  simulation  it  is  assumed  that 
the  external  reformer  is  operated  at  an  optimal  operating 
point.  Therefore  the  composition,  velocity  and  tempera¬ 
ture  of  the  anode  gas  is  considered  to  be  constant. 

•  The  air  at  the  cathode  inlet: 

For  the  cathode  air  the  composition  at  the  inlet  is  con¬ 
sidered  to  be  constant  (79%  N2  and  21%  O2).  But  the 
inlet  temperature  (Tc?in)  as  well  as  the  inlet  velocity  (wc,in) 
can  be  varied.  The  temperature  distribution  is  significantly 
affected  by  both  values.  With  an  increasing  velocity  or  a 
decreasing  temperature  at  the  air  inlet  the  fuel  cell  is  cooled 
more  intensively.  As  a  result  the  polarisation  resistances 
and  thus  the  voltage  losses  rise. 

In  contrast  to  the  characteristics  of  the  fuel  the  param¬ 
eters  at  the  cathode  inlet  can  be  changed  relatively  easily. 
The  velocity  is  specified  by  the  blower,  whereas  the  tem¬ 
perature  is  determined  by  a  heat  exchanger. 

•  Electrical  input  parameter: 


-0-  average  power  density 
— cell  voltage 
voltage  loss  by: 

— polarisation  resistance  of  02 
a  ohmic  resistance  at  the  cathode 

- transport  resistance  of  the  ions 

v  ohmic  resistance  at  the  anode 
— polarisation  resistance  of  H2 
-9-  polarisation  resistance  of  CO 


Fig.  10.  Resistance. 
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Fig.  11.  Operating  range  for  an  average  current  density  of  iei,avg  = 
250.0  mA  cm2. 

For  the  definition  of  an  operating  point  either  the  cell 
current  (galvanostatic)  or  the  cell  voltage  (potentiostatic) 
must  be  given. 

The  possible  operating  range  is  limited  by  material  prop¬ 
erties,  especially  the  maximum  thermal  stress.  For  this  sim¬ 
ulation  a  maximum  temperature  in  the  solid  of  Ts ,max,crit  = 
1250  K  is  considered.  Furthermore  the  average  solid  temper¬ 
ature  should  not  drop  below  rs>avgjCrit  =  1100K  or  the  cell 
voltage  decrease  to  a  point  where  the  nickel  catalyst,  which 
is  used  in  the  anode  electrode,  starts  to  react  electrochemi- 
cally.  At  a  temperature  of  Ts  =  1 123  K  a  critical  cell  voltage 
of  f/pc  =  0.6  V  should  be  considered. 

In  Fig.  11,  the  operating  range  is  displayed  for  a  current 
density  of  /ei,avg  =  250.0  mA  cm2  depending  on  the  cathode 
inlet  velocity  and  temperature.  The  boundaries  of  the  admis¬ 
sible  operating  range  given  by  the  average  and  maximum 
solid  temperature  are  also  drawn. 

Further  restrictions  for  the  operating  range  are  the  maxi¬ 
mum  temperature  gradient  in  the  solid  as  well  as  the  maxi¬ 
mum  blower  power  and  the  size  of  the  heat  exchanger.  Simi¬ 
lar  analysis  can  be  carried  out  for  other  current  densities  and 
thereby  an  optimal  operation  point  can  be  determined  within 
the  resulting  three  dimensional  space  at  which  a  maximum 
power  density  is  attained. 
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5.  Conclusions 

A  one-dimensional  stationary  model  of  a  SOFC  used  in  an 
Auxiliary  Power  Unit  has  been  presented.  At  the  cathode  of 
the  fuel  cell  the  reduction  of  oxygen  is  considered,  while  at 


the  anode  both  the  oxidation  of  hydrogen  and  the  oxidation 
of  carbon  monoxide  are  regarded.  For  the  polarisation  resis¬ 
tances  the  approach  of  Achenbach  [1]  is  used.  The  modelling 
covers  the  water-gas  shift  reaction  as  well  as  the  methane 
steam  reforming  reaction  in  the  anode  channel  of  the  fuel 
cell. 

The  commercial  mathematical  tool  FEMLAB  was  used 
to  implement  the  equations  and  to  calculate  the  results.  It 
provides  a  graphical  user  interface  for  the  implementation 
and  for  the  analysis  of  the  results.  Further  on  predefined 
application  modi  could  be  used  to  implement  the  model 
equations. 

A  suitable  operating  range  of  the  solid  oxide  fuel  cell  un¬ 
der  variation  of  the  inlet  temperature  and  inlet  velocity  of  the 
cathode  gas  as  well  as  the  anode  gas  could  be  predicted.  For 
a  current  density  of  /el,avg  =  250  mA  cm2  exemplary  simula¬ 
tion  results  are  shown.  The  co-flow  mode  is  to  be  preferred 
because  of  the  high  temperatures  within  the  active  cell  area 
for  the  counter-flow  mode.  The  operating  range  is  determined 
in  terms  of  temperature  and  amount  of  cathode  air  used.  Fur¬ 
ther  on  temperature  limitations  are  considered. 
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